Preprint typeset in JHEP style. - HYPER VERSION 



CERN-PH-TH/2011-082, LAPTH-010/11 



The Cosmic Linear Anisotropy Solving System 
(CLASS) II: Approximation schemes 



Diego Blas^*, Julien Lesgourgues'*'^'^, Thomas Tram 

""Institut de Theorie des Phenomenes Physiques, 
Ecole Polytechnique Federate de Lausanne, 
CH-1015, Lausanne, Switzerland. 

^ CERN, Theory Division, 
CH-1211 Geneva 23, Switzerland. 

" LAPTh (CNRS - Umversite de Savoie), BP 110, 
F-74941 Annecy-le-Vieux Cedex, France. 

Department of Physics and Astronomy, 
University of Aarhus, 
DK-8000 Aarhus C, Denmark. 



Abstract: Boltzmann codes are used extensively by several groups for constraining cos- 
mological parameters with Cosmic Microwave Background and Large Scale Structure data. 
This activity is computationally expensive, since a typical project requires from 10^ to 10^ 
Boltzmann code executions. The newly released code CLASS (Cosmic Linear Anisotropy 
Solving System) incorporates improved approximation schemes leading to a simultaneous 
gain in speed and precision. We describe here the three approximations used by CLASS for 
basic ACDM models, namely: a baryon-photon tight-coupling approximation which can be 
set to first order, second order or to a compromise between the two; an ultra-relativistic 
fluid approximation which had not been implemented in public distributions before; and 
finally a radiation streaming approximation taking reionisation into account. 
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1. Introduction 



Boltzmann codes have experienced considerable improvements in terms of precision and 
speed with respect to the pioneering COSMICS package |Q. In each new public code (CMBFAST 
||2[, CAMB Q, CMBEASY Q), several sophisticated optimisation methods and approximation 
schemes have been introduced. Efforts on this side keeps being justified for two reasons. 
On the one hand, we need to fit data with higher and higher precision. For instance, 
the analysis of Planck data requires much more accurate theoretical predictions than for 
WMAP On the other hand, a growing number of cosmologists are interested in fitting 
cosmological data with several extensions of the minimal cosmological model, in order to 
probe new physics. This requires running parameter extraction algorithms on computer 
clusters, often involving 10"^ or 10^ Boltzmann code executions (for each new model or new 
combination of data sets). Hence, any way to speed up Boltzmann codes without loosing 
precision is useful. 

In front of such needs, a new code, the Cosmic Linear Anisotropy Solving System 
(class) has just been released^. The goal of this project is not just to improve speed and 
precision, but also to provide a flexible and user-friendly code that can be easily generalized 
to non-minimal cosmological models. In this paper, we do not discuss flexibility issues and 
only concentrate on the improved approximation schemes used by CLASS, in the strict 
context of minimal ACDM cosmology. Extensions requiring extra approximations may be 
introduced and discussed case by case in the future. In a companion paper we already 
discuss the approximation specific to massive neutrinos and non-cold dark matter relics. A 
comparison between the power spectra obtained by CAMB and CLASS for the minimal ACDM 
model, as well as estimates of the relative speed of the two codes, is presented in |8|. 

The next three sections describe: a baryon-photon tight-coupling approximation which 
can be set to first order, second order or to a compromise between the two (Sec. |2|); an ultra- 
relativistic fluid approximation which had not been implemented in public distributions of 
Boltzmann codes before (Sec. |3|); and a radiation streaming approximation consistently 
including reionisation (Sec. Appendix ^ describes the stiff integrator which can be 
used by CLASS as an alternative to the Runge-Kutta integrator: without this integrator, it 
would have been essentially impossible to launch test runs with no approximation schemes, 
to evaluate the error induced by these schemes; moreover, this integrator gives better 
performances even in the presence of approximations, and is set to be the default integrator 
in CLASS. We summarize our full approximation landscape in Sec. ^. 

Note that throughout this paper, when discussing CAMB and CLASS, we refer to the 
versions available at the time of preparing this manuscript, i.e. the January 2011 version 
of CAMB and vl . 1 of CLASS. 

2. Tight Coupling Approximation (TCA) 

Before recombination, when the opacity aUeCFx is very large, the equations governing the 
tightly coupled baryon-photon fluid form a stiff system. Indeed, the opacity defines a 

^available at http://class-code.net 
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conformal time scale of interaction Tc = (aneor)"^ considerably smaller than that on which 
most of the modes actually evolve, namely th = a/ a' for super-Hubble scales and = 1/k 
for sub-Hubble scales (see Fig. |l|) . Standard integrators like Runge-Kutta algorithms would 
be very inefficient in solving such a system. This motivated Peebles & Yu to introduce a 
simplified system of differential equations valid in the regime of small Tc/th and Tc/t/; 
(tight coupling approximation or TCA) Q. The overall idea is that quantities which are 
vanishingly small in the limit Tc — )• are solved perturbatively in Tc, and these analytical 
expressions are used in the numerical code solving the remaining differential equations of 
the system. 



1000 f 1 , 1 1 1 , 6 
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Figure 1: (Left) Evolution characteristic conformal time scales in units of Mpc. Before recombi- 
nation, the baryon-photon interaction time scale Tc = {anecrT)~^ is much smaller than the Hubble 
time scale th = a/a'. For each mode and in a given range, it is also smaller than the acoustic 
oscillation time scale 1/k. (Right) Evolution of the product a scale factor times baryon adiabatic 
sound speed), with an arbitrary normalization of the scale factor, before and during hydrogen re- 
combination (in this model, the visibility function peaks at Tree = 278 Mpc). For r < 200 Mpc, 
helium recombination leads to a variation of this product by approximately 10%. 



The CLASS user can choose between two integrators for the system of linear perturba- 
tions. One of them (ndfl5, see Appendix^) is optimised for stiff equations, and shows 
good performances even in the tight coupling regime. However, by reducing drastically the 
number of equations to integrate, any TCA scheme leads to a speed up even in presence of 
such an integrator. 

Tight coupling equations have already been derived and improved by many authors 
after Peebles & Yu's seminal paper. For instance, the TCA formulas presented in Ma & 



Bertschinger |10| are derived to first order in Tc (omitting some polarisation terms which 
contribute to the photon shear at this order). Lewis et al. implemented the full first-order 
solution in CAMB also relaxing Ma &: Bertschinger 's assumption that Tc oc a^. Doran 
implemented some improved formulas in CMBEASY which are valid in the Newtonian gauge, 
and include a few contributions beyond order one [p] ]. 

As we will see below, the first-order TCA formulas provide poor approximations to 
the baryon-photon differential energy flux and to the photon shear at large times. This 
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is not much of a problem if one switches to the exact^ equations early enough. However, 
a better scheme would allow to switch off the TCA later, and to save a lot of integration 
time without loosing precision. 

Here, we will derive the full second-order TCA formulas in the synchronous gauge. 
While this work was in preparation, Cyr-Racine and Sigurdson published a paper on exactly 
the same topic [12|. We will show that the numerical results from our approach coincide 



with those of Another recent paper discussing the TCA beyond first order and its 
implementation in second order cosmological perturbation theory is ||l^ : this work actually 
proposes a systematic way to compute high-order corrections in a given model and at any 
order^. 

2.1 Full equations 



In the following we will adopt the notation of Ma & Bertschinger |1C]. The baryon pertur- 
bations will be characterized by the energy density contrast 6^ and the divergence of the 
fluid energy flux 6),. These quantities satisfy the equations^ 

K = - \h', (2.1a) 
= -m + c^fc^^b + ^e^fe, (2.1b) 



with % = aH = a' /a, where a is the scale factor, and we have defined Q^b = 6^ — 9b and 
R = The quantity G^b represents the divergence of the energy flux of the photons 
(^-y) in the frame comoving with the baryons. Its time-derivative Q'^^ is often referred to 
as the "baryon-photon slip". The fields h and r] represent the metric perturbations. As 
explained in |l^ , the baryon pressure perturbation can be safely neglected in the continuity 
equation, but its Laplacian should be kept in the Euler equation (giving raise to the term 
c^k'^Sb) since it affects the evolution of very small wavelengths, smaller than the baryon 
Jeans length. In [l^ and many other references, the baryon sound speed is identified to 
the adiabatic sound speed 

,_kBTb (2.2) 



/i \ 3 din a 
where the evolution of the proton temperature is given by 

Ti = -2nn + n) , (2.3) 

nieTc 

and /i is the mean molecular weight. This approximation has been proved to be inaccurate 
in refs. [|15|, [iql, but the difference is only important for computing the matter power 



^It should be intended that throughout this section, the word "exact" is intended in the sense "without 
using the TCA". Obviously, our equations are never exact since they rely on a number of common and 
well-justified approximations, e.g.: linear perturbations, pressureless CDM, approximate expression for 
photon-baryon coupling, baryon pressure neglected in several equations, etc. 

^As pointed out by the author of [lij, this method can be implemented numerically provided that the 
code computing the evolution of thermodynamical variables outputs fully continuous and derivable functions 
of time, which is precisely the case with the CLASS version of RECFAST. 

*We will work in Fourier space and stick to the synchronous gauge in conformal time, which we will 
denote by r. We use the prime to denote derivative with respect to conformal time. 
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spectrum for k ^ 10 /iMpc^^. The current version of CLASS (vl . l) and CAMB (from January 
2011) still relies on the = approximation, while future versions of both codes are likely 
to switch to the actual sound speed calculation, as in the CAMB_source^ code. This issue 
is irrelevant for the results of this paper, which do not involve very large wavelengths. All 
numerical results below have been obtained using instead of c^, but our formulas can 
adequately describe the large k range, provided that a correct sound speed calculation is 
performed. 

Like the adiabatic sound speed, the actual baryon sound speed is expected to decrease 
approximately as o~^, except during helium and hydrogen recombination (see Fig. [l|). 

The characterization of the photon distribution requires the determination of its differ- 
ent multipoles 6.y, 9^, a-y and F^i for / > 3. These satisfy the recursive Boltzmann equations 
(Eqs. (63) in P): 

6'^ = -^e^ - (2.4a) 
e'^ = km6,-a,)-^, (2.4b) 
K = ^^7 - f ^^73 + Mh' + 67?') - ^CT^ + ^ (G,o + G,2) , (2.4c) 
= 2ITT [^^7(/-i) - + l)^7(m)] - :^^7^ ^>3 (2.4d) 
^'-ri = Wn [^'^7{«-i) - + l)G'7{'+i)] 

+^ [-G^i + i {F^2 + G^o + G^2) [Sio + (2.4e) 

where we remind the reader that Ky2 = 2a-y. Finally, the previous hierarchical equations 



are truncated at some I = lma,x following equation (65) in |1C] 



T 

2.2 TCA equations 



From Eqs. (2.1) and (p.4|) one sees that the different time scales in the problem are Tc, 



k ^ and the time scale of cosmological evolution H. As we will show explicitly in the 

)'tca tc: 
-yb ' "7 



next section, it is possible to find a solution (Qy^^, cr*'^^) for the baryon-photon slip and 



the photon shear in terms of 6h, 6^, Oh and 6^, which is valid to any desired order in the 
small parameter Tc. The knowledge of this solution helps to reduce the full system of 



equations (2T) and ( |2.4[) to just four of them for the low multipoles of the distributions. 



More concretely, one may use Eq. ( p. la ) for 6h' , and Eq. ( p. 4a ) to determine 6y'; the energy 



fluxes are characterized by the linear combination of Eqs. ( 2.1b| ) and ( p. 4b ) in which the 
coupling term vanishes: 

eb' + R9^' = -n9b + clk^Sb + Rk^(^^6^-a'^'''^ , (2.6) 

and finally O-y' — 6b = &jb^- As desired, this scheme allows us to get rid of any coefficient^ 
in T~^. For practical reasons it is customary to combine linearly the last two equations 



^http : / / camb. info/ sources/ 

®The scale Tc appears now as a small perturbation. 
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in order to eliminate QJ and get an expression for only; then can be found from 



equation (|2.6| ). In summary, once Ql^'^ and o^^^ are known the goal is to solve the closed 
system formed by the four equations ( |2.1a ), ( 2.4a| ), and 



Udi, - ctk'^5b - k'^R i T<57 - <^ ) + RQ 



(1 + R) 



)'tca 
76 



(2.7a) 
(2.7b) 



2.3 Perturbative expansion 

The aim of this section is to find expressions for 0^^^ and valid at the ii}^ order in Tc 
(the zero order is trivial: both species behave as a single perfect fluid, so that and all 
multipoles of the photons beyond 5^ and 9^ vanish; the first order can be found in [0]). 
We first multiply equations ( 2.4b|) and ( 2.1b| ) by Tc- 

9^' - [j^i - ^^7) + = , (2.8a) 
Tc [-Ob' - ndb + c^k'^db] + RQ^b = . (2.8b) 
To obtain a differential equation for Q^b we can combine the above two equations into 



Q^b — 'H9b + k^ { (?s6b — + o"^ 



+ (1 + R)@^b = . 



(2.9) 



This equation involves the photon shear, given by the following equation (see ( ^.4c| )): 
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4w 



-9^ + -h' + 8r/' - lOfj^' - 3A;F^3 



+ ^^(^70 + ^72). (2.10) 



Until now, all these equations are exact. In the limit of interest, both of them can be 
schematically written as'' 

ey{t)' + y{t)/f{t) + eg{t) = Q, (2.11) 

where e is a small parameter. In our case, e can be chosen to be fc, the opacity at an 
arbitrary time around which the expansion is performed®. The perturbative solution is 
given by 

y(t) = ^e"y„(t), yi = -fg, = -M . (2.12) 

n=l 

Notice that for functions f{t) and g{t) with smooth time variations on the scale Tc, the 
previous is a perfectly well defined solution. Finally, the most general solution is found by 
adding to the previous particular solutions the solution of the homogeneous equation 



ey{ty + y{t)/f{t) = 0, y = Ce-^'^^r'^' . 



(2.13) 



^This is immediate for (^.£|). For ( |2.10| ) if follows from (2.4e), and we will explicitly verify it shortly. 
*In fact, the small dimensionless parameters will be fck and TcH, with T-L evaluated around the same 
time as tn. 
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Note that in our case / is always positive, which is enough to make this part of the 
solution suppressed very fast. Hence, the relevant solution is given by Eq. (2.12), which, 
after absorbing the small parameter in the function f = ef, reads: 



y{t) = ^yn{t), yi = -fg, yn+i = -fy'n 



(2.14) 



n=l 



In terms of these functions, the powers of / and its derivatives represent the different orders 
of the approximation. 

2.4 Second-order approximation 

Using the previous expansion in Eq. ( |2.g| ), the baryon-photon relative velocity reads at 
order two: 

= /e (-50 + lege + feg'e) + 0{f!) , (2.15) 

with 



fe = rr^, ge = -'Heb + k'{ci5b-i6^ + a^] . (2.16) 



1 + i? ' 

We still need to differentiate this equation in order to get a similar approximation for the 
slip: 



el 



76 



Ik 

fe 



Qyb + fe (-g'e + fe9e + 2/45e + /e5e ) + ^(r^) . (2.17) 



There are several ways to organise and simplify the final result. In order to write the 
first-order term in the same form as in the rest of the literature, we need to express g^ in 
a very peculiar way: 

g'e = -ne', - n'Ob + (^cf 4 + 0% - -^5'^ + a'^ 

{Ucl + cf )(5b + c% - -<5; + a'A + —Q^b (2.18) 



2^eV ^fe + 

a 



n 



7 ^ ^^7 



(2 + R)n 



G 



7b 



Here we defined = {7ic^ + cf): this quantity would vanish if the approximation oc 
were valid at all times. In the second line, we used Eq. ( 2.1b| ), while in the third line we 
used ( 2.4b| ): so these expressions for are all exact. 

The first-order approximation for 0^^^ is obtained by replacing the first occurrence of 
g'e in Eq. ( 2.17] ) with the last expression of (2.18), in which we neglect the terms 
and {2%(T^ + a'^) which represent contributions of higher order. The final result is: 



61 
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< 2% 

Tc l + R 



e^b - fe ^--Ob + [-jS^ + c% + c% - + 0{fl). 

(2.19) 

Note that in the previous expression we used the exact relation R' = —TiR. 

For the second-order expression for the slip, we go back to Eq. ( |2.17 ). We replace the 
first occurrence of g'^ by the full expression ( |2.18|) , assuming that Q'^^ and {2%a^ + a'^) 
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have been replaced by their first-order approximation. Finally, we replace g@, q'q and Qq in 
the last three terms by their zeroth-order approximation. The final result can be written 
in a compact form: 



GL = (l-2?^/e) U 



2% 



T, l + R 



©76 - /e 



n 



1./ 



-fek'' {2na^ + a;) + /e [f^ge + 2f'^g'^ + ~feg'^\ + 0{fl). (2.20) 

This formula requires an expression for the shear valid at order one. However, to solve 
equations ( 2.1b , 2.4b ) consistently to second order, we need the expression for the shear at 
the corresponding order. This can be achieved using ( |2.1ClD . To solve this equation let us 
first note that the polarisation multipoles i = 0, 2 obey (cf. ( |2.4g| )) 



70 



kG^i + 



-G^o + C7 + - (G^o + G^2) 



, k I 
= g (2G^i — SG^s) + 



-G^2 + (2o"^ + G-yo + G^2) 



(2.21) 



from which we see that, at first order in Tc, G 



72 



G 



70 



cr^. Thus, it is consistent to 



consider these multipoles as 0{tc), and write the second order solution to ( 2.10 ) following 
(I2II as 



-e^ + -h' + 877' 



IOtp d /t, 



dr V 9 



8 4 



+ ^[G,o + G^2]) +0(Tc3), 



(2.22) 



where we used the fact that F^s = O(tc^). Indeed, the high photon multipoles obey at 
leading order (cf. ( |2.4d )) 



iTck 

21 + 1 



The final step involves an evaluation of ( ^.211 ). Again, from ( 2.4e ), one finds that 

G^i ~ G^3 ~ O(tc^). 



This allows us to find the perturbative solution to ( 2.21 ) 

G7O = — - ^TcCT^ -M^I^Tc Lt^2 

From the previous expression and (|]2|) we find 



which implies 



8t^ 
45 



a, = ^(20^ + /i' + 6r/')+O(f2), 



20; + /i" + 6r?") + ^ (20^ + /i' + 6r/') + 0(f2). 



(2.23) 



(2.24) 



(2.25) 
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Finally, the shear at second order is found from ( 2.22 ) to be 

llr, 



87^ 

45 



[26^ + h' + 6r/') 



11^ 
6 



6 



(2^; + h" + 67/") 



+ 0(f3). (2.26) 



The last missing items are the zero-order expressions for qq, and appearing in 
equation ( |2.2C1| ). Noticing that equation ( p. 7a ) implies 



1 



l + R 



-nOb + k^ci6h + k'R-6j ) + 0{fc), 



RH „ ^ r-/ 



+ 0(tc 



we can write these last terms as 



9e 



9'e 



9'^ 



0{fc), 



(2.28a) 
(2.28b) 



-ne'^ - 2n'ei - n"eb 



+ k' 



{cl)"6b-2{ci) 



2\i 



+ 0(fc). (2.28c) 



The derivation given in [|T^] follows different steps, but since it is still a second-order TCA, 
the results should be identical under the approximation oc used in [12|, at least up to 
terms of order three or higher. For the shear, our expressions are indeed exactly identical. 
For the slip, there are so many ways to write the result and so many terms involved that the 
comparison is not trivial. However, by coding the two formulas in CLASS and comparing 
the evolution of Q^b in the two cases, we found that the two expressions agree very well, 
since numerically the difference appears to be at most of order 0{f^). 

2.5 Implementation of various schemes in CLASS 

We implemented various TCA schemes in CLASS, which can be chosen by setting the 
flag tight_coupliiig_approxiniation to different values. We always start integrating the 
wavenumbers very deep inside the tightly-coupled regime. Hence, in contrast to [12|, we do 
not need to include terms in in the initial conditions. We always use the set of equations 
(|2li]) , (|1D, (|7|), (|27b| ) with different expressions for the slip Q^^^ and the shear a^^: 



oc a 



1. first-order expressions ( p.l9| , 2.24), with the approximations Tc oc and 
used e.g. in ref. [10]. This corresponds to the scheme used by CLASS when the label 
tight_coupling_approximation is set to f irst_order_MB. 

2. first-order expressions ( |2.19 , 2.24| ) with the only approximation oc a~^, like in 
CAMB. This scheme is used by CLASS when the same flag is set to f irst_order_CAMB. 

3. exact first-order expressions (|1|, |]2|) when the flag is set to f irst_order_CLASS. 

4. second-order expressions from [|l2| when the flag is set to second_order_CRS. 
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5. second-order expressions from Eqs.( ^.20 , 2.26D for the flag second_order_CLASS. 



6. finally, second-order expression for the shear, but a reduced expression for the slip, 
involving only the leading order-two terms: 



76 

2 



fek' 



(2.29) 



This option is taken when the flag is set to compromise_CLASS, and is chosen to be 
the default option in CLASS. Notice that the last scheme does not have a term h" , 
which is advantageous from the computational point of view^. 



To justify the compromise scheme, notice that it encapsulates the leading order in ( 2.2^ ) 



for subhorizon modes (note that has an extra momentum dependence with respect 
to the other perturbations, and each time derivative adds one more power of k in this 
regime). Thus, from Eq. ( p.20| ) we learn that Eq. (2.29) implements the leading second 



order correction for the modes with a big comoving momentum. It is precisely for these 
modes that the first order approximation fails first, which explains the success of the 
compromise scheme. To check that the approximation is indeed correct, we implemented 
this scheme for several modes k assuming ACDM. As illustrated in the next subsection, 
this scheme is nearly as good as the full second-order one, being at the same time much 
more compact and requiring many less floating point operations. 

2.6 Comparison at the level of perturbations 

In Fig. ^, we compare these different approximations for a flxed wave number (chosen to be 
k = lO^^Mpc^^), and in Fig. |3| for a flxed conformal time r (chosen to be the time when 
the TCA is switched off in the previous example). 

Scrutinising first the various first-order schemes, we see a significant difference at late 
time between the first two (MB and CAMB), showing that TcCC is a bad approximation. 
However, there is no sizable difference between the second one (CAMB) and the third one 
(class) in which the approximation oc is relaxed. We reach the same conclusion 
when comparing second-order schemes with or without the same approximation. This is 
not a surprise, since we are only considering scales larger than the baryonic Jeans length 
at any time. When studying very small wavelengths, the CLASS user is free to choose one 
of the TCA schemes where the full evolution of is automatically taken into account 
(namely, f irst_order_CLASS or second_order_CLASS), but as mentioned before, in this 
limit, an accurate sound speed computation should also be implemented in order to relax 
the (?g = approximation. 

Figs. § and ^ show that all first-order schemes provide poor approximations for the 
slip and the shear near the end of the tightly coupled regime, and hence, inaccurate initial 



®In order to compute h" one should use one more Einstein equation that in the standard case, and 
compute the 5TI component of the stress-energy tensor, i.e. the pressure perturbation for all species. 
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Figure 2: Evolution of 6-,^ (left) and cr^ (right) for the mode k — lO^^Mpc^^, using the various 
TCA schemes listed in Sec. In each case, the quantities are represented as points when the TCA 
is switched on, and as continuous lines of the same color when exact equations take over. The TCA 
is switched off at r = 1 Mpc in the reference case, and at r = 194 Mpc in all other cases. We show 
a single set of points for cases which are indistinguishable by eye, namely: f irst_order_CAMB and 
f irst_order_CLASS; and also, second_order_CRS and second_order_CLASS. The default scheme 
compromise_CLASS is also hardly distinguishable from the second_order_CLASS. 
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Figure 3: Comparison of <d^i,{k) (left) and cr^{k) (right) at the time when the TCA is switched 
off (r — 194 Mpc in this case), for the different schemes listed in Sec. 2.5. All cases are compared 
to the full second-order scheme second_order_CLASS. For Q^b^ we show a single curve for the 
cases f irst_order_CAMB and f irst_order_CLASS, since they are indistinguishable by eye. For 
cr-y, all schemes using the first-order shear expression are indistinguishable; this is also true for all 
cases using the second-order shear expression (despite the fact that 9~f and metric perturbations are 
slightly different in each individual case). 



conditions at the time at which the full equations are turned on. As expected, second-order 
schemes work much better. We find a very good agreement between secoiid_order_CRS 
and second_order_CLASS: this validates both the results of |12] and our results. The 
residual difference is likely due to the fact that our two independent derivations lead to 
expressions in which some higher-order terms (of order 0(f^)) appear in different ways. 

These two schemes also agree to a very good extent with compromise_CLASS, which 
is much more straightforward to code, and computes the baryon-photon slip with approx- 
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imately ten times less operations. In particular, with this scheme, the code does not even 
need to compute the quantities h" , t" and 7i", which are not so obvious to obtain with few 
operations and good accuracy. Hence, this method is set to be the default TCA in CLASS. 

2.7 Comparison at the level of temperature/polarisation spectra 

First, let us specify which precision parameters in CLASS govern the evolution of perturba- 
tions in the early universe and the TCA switching time: 

• two parameters define the time at which initial conditions are imposed during the 
tightly-coupled stage. Each wave-number starts being integrated (with one of the 
TCA schemes) as soon as one of the two conditions 

{tc/th) > start_small_k_at_tau_c_over_tau_h 

or 

(r/z/rfc) > start_large_k_at_tau_h_over_tau_k 

is fulfilled. The second condition means that at initial time, wavelengths should be 
sufficiently far outside the Hubble scale; the first condition, which over-seeds the 
second one for the smallest wave numbers, means that the initial time should not be 
too close to recombination. 

• two parameters define the time at which the TCA is turned off for each wave number. 
This happens when one of the two conditions 

{tc/th) > tight_coupling_trigger_tau_c_over_taiL_h 

or 

i'Tc/'Tk) ^ tight_coupling_trigger_tau_c_over_tavL_k 

is fulfilled. CLASS imposes that the TCA switching time should always be chosen after 
the initial time, which means that the four parameters above should satisfy simple 
inequalities. 

• one parameter defines the time (common to all wave numbers) at which the source 
functions (leading to the computation of temperature and polarisation C;'s) start 
being sampled and stored. This happens when the condition 

(tc/th) = start_sources_at_tau_c_over_tau_h 

is satisfied. This time can eventually be chosen during the tight-coupling regime for 
the smallest wave numbers. 

In order to show the impact of these parameters, we take the set of precision parameters 
defined in the file cl_permille .pre of the CLASS public distribution, which corresponds 
to an accuracy of at least 0.1% on each temperature and polarisation Ci, and uses the 
coniproniise_CLASS scheme. We then vary the two trigger parameters mentioned above, as 
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described in Table ||. The first setting, called no-tca, corresponds to switching off the tight 
coupling approximation immediately after setting the initial conditions, so that no TCA is 
ever used. This leads to reference temperature/polarisation spectra with respect to which 
all the other results of this section are compared. The settings called teal, tca2 and tcaS 
introduce from 0.02% to 0.08% of error with respect to the no-tca case, as illustrated in 
Fig. I 

In Fig. |5|, we stick to the precision setting tca3 and compare the different TCA schemes. 
The f irst_order_CAMB and f irst_order_CLASS results are indistighuishable, confirming 
the fact that the approximation oc a^^ is sufficient in practice. Our second-order results 



and those derived from [12| are also in perfect agreement. As expected, the results from the 
compromise_CLASS scheme are essentially as good as the full second-order results, while 
the first-order results are roughly ten times less accurate. We also show on this plot the 
error produced by the f irst_order_CLASS scheme with teal precision settings, which is 
similar to that produced by the compromise_CLASS scheme with tca3 precision settings. 
Hence, in order to estimate the usefulness of going beyond the first order TCA, we can 
compare the performances of the code in these two cases. 

In Table ^, we compare running times in the no-tca case and in the previous two cases, 
using either the Runge-Kutta or ndf 15 integrator. The timings displayed here correspond 
to the number of seconds spent by our computer in the perturbation module of CLASS, in 
a non-parallel execution. The ndf 15 integrator is always better, by a huge amount in the 
no-tca case, or by 20 to 30% in the other cases. Being unaffected by the issue of integrating 
a stiff system, the ndf 15 integrator is not very sensitive to the choice of TCA scheme, with 
only a 3% speed up when using the compromise scheme instead of first-order schemes. The 
Runge-Kutta integrator is more sensitive, with a 9% speed-up for the compromise scheme. 

We conclude that CLASS benefits much more from the implementation of our stiff 
integrator than from going beyond the first-order TCA. For some particular models, the 
user may wish to stick to the Runge-Kutta integrator, in which case the compromise_CLASS 
scheme leads to a sizable speed up. 





no-tca 


teal 


tea2 


tea3 


tight _coupling_trigger_tau_c_over_tau_h 
tight _coupling_trigger_tau_c_over_tau_k 


4.1 • 10"^ 
6.1 • IQ-^ 


7-10-3 
3 • 10-2 


8 • 10-3 
5 • 10-2 


9 • 10-3 
8 • 10-2 



Table 1: Four settings for the precision parameters governing the time at which the TCA is 
switched off. 



3. Ultra-relativistic Fluid Approximation (UFA) 

All massless neutrinos and ultra-relativistic relics can be treated as a single species, labeled 
as "ur" in CLASS. The code assumes that these species are fully decoupled. Hence, they just 
free-stream within a given gravitational potential, and can be followed with the collisionless 
Boltzmann equation expanded in harmonic space and integrated over momentum [^] . The 
solution can be formally written as the sum of spherical Bessel functions j; (kr) (exhibiting 
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Figure 4: Impact of precision parameters governing the TCA switching time, for temperature 
(left) and E-polarisation (right). We show the power spectrum of the three settings teal, tca2 and 
tcaS compared to the reference spectrum no-tca (see Table |^ for precision parameter values). 



precision 


no-tca 


teal 


tca3 


TCA scheme 


(irrelevant) 


f irst_order_class 


compromise_class 


rk 


1069s 


19.4s 


17.8s 


ndf 15 


16s 


14.9s 


14.6s 



Table 2: Execution time of the perturbation module, in seconds, with several TCA settings and 
with the two integrators (Runge-Kutta and stiff integrator ndf 15). The last two columns lead to 
roughly the same level of accuracy. 



damped oscillations for r > l/k), plus a particular solution of the inhomogeneous equations 
sourced by metric fluctuations. The fact that one part of the solution has an analytical 
expression cannot be directly implemented in the code, because the total perturbations 
(^uri ^ur and cJur back-react on the metric perturbations through Einstein equations, and 
affect the source terms in the Boltzmann equation. We will use this decomposition only as 
a guideline for deriving accurate approximation schemes. 

3.1 Truncation of the Boltzmann hierarchy 

Since the ur species couple only gravitationally to other species, we are only interested 
in tracking Jur; ^ur and cTur- Higher multipoles must still be included since they couple 
with the lower ones, but in all efficient Boltzmann codes, the hierarchy is truncated at 
some low multipole value /max- CMBFAST, CAMB, CMBEASY and CLASS all use the truncation 
scheme proposed in Ma & Bertschinger (Eq. (51) of [^]) which is designed to minimize 
artificial reflection of power from Zmax back to lower multipoles. Still, this truncation is 
not perfect, and a significant amount of unphysical reflection cannot be avoided for times 
beyond r = /max/^- This implies that in order to compute an accurate CMB spectrum, /max 
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Figure 5: Impact of TCA schemes on temperature (left) and E-polarisation (right). For the 
various TCA schemes discussed in Sec. ^.Sl , wc show the power spectrum with accuracy settings 
tcaS compared to the reference spectrum no-tca (see Table [l] for precision parameter values). We 
show a single curve for cases which are indistinguishable by eye, namely: f irst_order_CAMB and 
f irst_order_CLASS; and also, second_order_CRS and second_order_CLASS. The default scheme 
compromise_CLASS is also hardly distinguishable from the second_order_CLASS. The faint line 
shows for comparison the first-order results with accuracy teal: the error is then comparable to 
compromise_CLASS with tcaS. 



should be at least of the order of 30. The computation of the matter power spectrum P{k) 
on small scales, up to some wavenumber fcmax; is more problematic: one should further 
increase /max proportionally to /cmax in order to get converging results. 

3.2 Sub-Hubble fluid approximation 

The Ultra-relativistic Fluid Approximation (UFA) implemented in CLASS is based on the 
idea that for a given wavenumber, /max should not necessarily be fixed throughout the 
whole time evolution. The code considers two regimes: wavelengths larger or compara- 
ble to the Hubble radius, and wavelengths much smaller than the Hubble radius. The 
transition between the two regime occurs for each k when the product kr (equal to 
r/rfc and coinciding with Tn/rk during radiation domination) reaches some threshold 
value that we call here (A;T)ufa. Typically, (A;r)ufa is chosen in the range from 10 to 
50, depending on the required precision. The full name of this parameter in the code 
is ur_f luid_trigger_tau_over_tau_k. In the first regime kr < (fcr)ufa, the Boltzmann 
hierarchy can be truncated at some /max which can be chosen to be rather small: it is 
enough to take to /max ~ (fc''")ufa! since multipoles with / > kr are negligible (according to 
the spherical Bessel function approximation). In the second regime kr > (A;r)ufa, the code 
still follows the three variables 6^1, O-ar and a-ar, which are sourced by metric perturbations. 
But multipoles in the range 2 < / <C fcr are suppressed, leading to an effective decoupling 
between the first three multipoles and the highest ones. Hence it is natural to lower /max 
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down to two in this regime. Ultra-relativistic neutrinos are then described by a reduced 
system of equations for 6ut, Our and cJur, i-e. by fluid equations (of course, this fluid is not 
assumed to be perfect, since it has anisotropic pressure). In summary, the UFA approxi- 
mation consists in lowering /^ax from a value close to (fcr)ufa down to /max = 2 deep inside 
the Hubble radius, at the time when kr = (A;r)ufa. Such a scheme offers many advantages: 

1. When computing the matter power spectrum, the number of ur equations to integrate 
before the approximation is switched on does not need to be scaled linearly with the 
highest wave number /cmax- 

2. The number of ur equations reduces to (/max + 1) = 3 in the whole region of {k,T) 
space fulfilling the condition kr > (fcT)ufa; this is precisely the region in which the 
computation would be time-consuming, since ur perturbations oscillate inside the 
Hubble radius. 

3. The UFA completely avoids the issue of power reflecting at some large /max; which 
would otherwise affect the evolution of low multipoles periodically due to some spu- 
rious wave travelling back and forth between / = /max and / = 0. (This behaviour is 
clearly seen in Fig. 

The fluid approximation could in principle be used until present time, but the code allows a 
more aggressive approximation, the Radiation Streaming Approximation, to take over from 
the UFA after photon decoupling. This new approximation is discussed in the next section. 
Hence, the UFA is essentially a way to save computing time during radiation domination 
and at the beginning of matter domination. 

3.3 Fluid equations 

We need to find a closed system for the evolution of , G-ar and , valid deep inside the 
Hubble radius. The full system of equations in the synchronous gauge (Eq. (49) in |]lO|) 
reads: 

•^ur' = -l^ur - (3.1a) 

Our' = k^S^ - a^) , (3.1b) 

2a^r' = ^^ur - ffci^urS + ^{h' + Gf]'), (3.1c) 
Krl=2m [^^nr^-D " (^ + l)Fur(l+l)] ■ (3.1d) 

In Appendix we use the formal solution of these equations in order to derive an exact 
integral relation between a^r'i Cur, ^ur and metric perturbations. We then find an approx- 
imate but more practical form of this relation valid inside the Hubble radius, at leading 
order in an expansion in metric perturbation derivatives {h^'^^/k"~^, r]^"'^ /k"'~^) and in 
powers of (fer)~^: 

3 2 1 

(^UT = CTur + ^^ur + ^/i' • (3.2) 

TOO 

Since metric perturbation only evolve over a Hubble time scale inside the Hubble radius, we 
expect this expansion to converge, and we will see below that the above relation is indeed 
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accurate enough for our purpose. In the default version of CLASS, this equation is used for 
closing the system of equations when the UFA is switched on. This method corresponds 
to the setting uf a_method = uf a_class in the code's precision parameter structure. 

3.4 Alternative schemes 

Some nearly equivalent schemes can be justified in slightly different ways. Truncating the 
Boltzmann equations at Imax = 2 with the usual truncation scheme of Ma & Bertschinger 
gives: 

(tJ = --(Tur + + l{h' + 6ri') . (3.3) 
T 6 6 

This truncation scheme is based entirely on the assumption that F^j- i{k,T) oc ji{kT). So, 



the reason for the difference between eq. (3.2) and (3.3) is that (|3.2D is based on the full 



formal solution, including the leading order contribution to the part sourced by the metric. 



while (3^) is based only on the solution of the homogeneous equation. Equation ( |3.3| ) is 
used when the user switches to ufa_niethod = ufa_mb, and amounts to adding an extra 
term in r}' . 



Finally, in a more general context, Hu |13] introduced a set of equations modeling a 



cosmological viscous fluid, and suggested that this fluid could approximate the evolution 
of free-streaming neutrinos with the parameter choice (w, c^, c^jg) = (1/3,1/3,1/3). In 
this limit, Hu's fluid equations are identical to our UFA equations except for the shear 
derivative: ^ 

^ur' = -3-a,r + l^ur + \{h' + 6r/') . (3.4) 
a 6 6 

The coefficient —3^ reduces to — 3/r deep inside the radiation dominated regime, but 
becomes different around the time of equality. This equation is used when the user switches 



to uf a_method = uf a_hu. Below we will compare the performances of equations (3^, p.3| . 



3.4) and show that the first one is slightly more precise (as expected from the rigorous 
mathematical proof of Appendix P). 

Finally, when the user selects uf a_niethod = ufa_none, no UFA scheme is employed, 
and the truncation multipole /max remains the same throughout the evolution. 

3.5 Comparison at the level of perturbations 

In Fig. |, we compare the evolution of Jur and for a given mode, obtained either by 
solving the full Boltzmann equation up to a very high /max ~ 3000, or with /max = 46 
with/without the default UFA. In absence of approximation, one can see some spurious 
evolution appearing periodically (here, around kr ~ 100 and kr ~ 200): this corresponds 
to the propagation of an unphysical wave between the multipole boundaries / = /max and 
/ = 0. Using the default UFA scheme uf a_class, we reproduce accurately the phase, the 
amplitude and, to a lesser extent, the zero-point of the oscillations. For the clarity of the 
figure, we do not show the results from alternative approximation schemes. We checked 
that the uf a_mb scheme also reproduces the correct phase and amplitude, but introduces 
a larger error in the zero-point of oscillations. Finally, the uf a_hu scheme reproduces the 
phase, but not the correct amplitude of the oscillations. 
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Figure 6: Evolution of (5ur (left) and a^r (right) for the mode k — over the range 

10 < fcr < 200, i.e. well inside the Hubble radius and throughout the radiation dominated stage. In 
the "exact" case, the full Boltzmann equation for ur is truncated at ^max ^ 3000, with no impact of 
the truncation on the result. In the "no approx." case, the truncation is performed at Zmax = 46. In 
the "CLASS" case, we use the UFA scheme uf a_class and set Zmax — 46 as long as kr < (fcr)ufa — 50, 
or Zmax = 2 afterward. 



3.6 Comparison at the level of CMB /matter power spectrum 



UFA scheme (fixed /max and (fer)ufa) 


uf a_none 


uf a_class 


rk 

ndf 15 


29.7s 
16.7s 


27.0s 
15.2s 



Table 3: Execution time of the perturbation module in seconds, with the precision parameters of 
the file cl_perinille .pre, plus Imax — (fc''')ufa — 18. Using an ultra-relativistic fluid approximation 
leads simultaneously to a 10% faster execution and to more accurate results. 

In Fig. 0, we show the impact of the UFA on the CMB and matter power spectrum. We 
take the precision parameters of the file cl_permille .pre, and play with the values of /max 
and (fcT)ufa, called l_max_ur and ur_f luid_trigger_tau_over_tau_k in the code. We 
first compute some reference spectra with /max = 3000 (to remove any truncation effect), 
and such a large value of (A;r)ufa that the UFA is never used. All other results from this 
section are compared to these spectra. We then fix both Zmax and (A;r)ufa to 18 and vary 
only the ur_f luid_approximation parameter. We show the error induced by each UFA 
scheme for the temperature and matter power spectrum in Fig. ^ (results for temperature 
and polarisation are very similar). The results from the ufa_none case are very unstable 
and depend a lot on the choice of /max and (A;r)ufa) values. With the present choice, they 
correspond to a twice larger error in the CMB spectra than in any UFA scheme; for slightly 
different choices they would also induce a larger error in the matter power spectrum. The 
three UFA schemes, which do not have such instabilities, are nearly as good as each other 
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Figure 7: Power spectra (for temperature CMB anisotropies and for matter) obtained with 
various implementations of the UFA, compared to those obtained in a reference run. In the four 
cases displayed here, we use ^max = 18 as long as fcr < (fcT)ufa — 18. The curves labeled "MB", 
"HU", "CLASS" correspond to the three possible implementation of the UFA discussed in the text; 
the last curve uses no UFA and a standard truncation scheme at Zmax = 18 until the time at which 
the next approximation RSA takes over (see Sec. ||): in this case the code is at the same time a bit 
less precise and 10% slower. We do not show the results for the polarisation spectrum Cf'^, that 
look very similar to those for C^^ . 



for CMB spectra, while for the matter power spectrum the uf a_class scheme is one order 
of magnitude better. Table ^ shows that the UFA approximation allows for a 10% speed 
up, while being more accurate for a fixed /max- 



4. Radiation Streaming Approximation (RSA) 



After their respective decoupling time, the photons gradually free-stream like neutrinos 
(except around the time of reionisation at which their coupling to baryons is enhanced). 
In principle, it would be possible to look for a fluid approximation for photons, like we 
did for neutrinos in the previous section. However we can go further than that, since 
during this period the universe is dominated by matter and eventually A/Dark Energy: in 
this case photons and massless neutrinos almost behave like test-particles in an external 
gravitational field, and we do not need to catch their evolution with high accuracy (which 
was not the case for the ur species during radiation domination). 

Like for massless neutrinos, in all efficient Boltzmann codes, the Boltzmann equation 
for photons is truncated at some low multipole value /max using the truncation scheme 
proposed in Ma & Bertschinger (Eqs. (65) in [10|). If /max is not large enough, the spurious 
reflection of power induced by the truncation propagates to the final results, because radia- 
tion perturbations still play a small role during the free-streaming regime. More precisely: 
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• the photon density fluctuation, velocity and shear perturbations appear in Einstein 
equations; 

• the photon density fluctuation and shear appear in the temperature/polarisation 
source functions; 

• the photon velocity appears in the evolution equations of baryons (since the baryon- 
photon coupling is not negligible during reionisation). 

In order to avoid propagating such an error, one can either increase /max, or find a way to 
infer the photon density, velocity and shear from some analytic Radiation Streaming Ap- 
proximation (RSA), in which case the integration of Boltzmann equations can be stopped 
soon after photon decoupling. Mathematically, this analytic approximation should coincide 
with the particular non-oscillatory solution of the inhomogeneous Boltzmann equations. 
Once the damped oscillations accounted for by spherical Bessel functions becomes negli- 
gible (i.e. when kr ^ I), the analytic approximation will coincide with the true solution. 
Before that time, it will provide a correct approximation to the true quantities averaged 
over a few oscillations. In summary, the RSA has a double goal: to avoid unphysical oscil- 
lations created by the Boltzmann truncation, and to avoid wasting time in integrating the 
Boltzmann equations over many such oscillations. 

The RSA does not need to be accurate at very late time (end of matter domination, 
A/Dark Energy domination), since by then radiation fluctuations are always completely 
negligible with respect to matter fluctuations. However, it should be reasonably accurate 
soon after photon decoupling, i.e. during matter domination, when the energy density in 
radiation, O,. = + ^lur, is smaller than one but not yet much smaller. The advantage 
of a better approximation is two-fold. First, it can be switched on earlier. Second, before 
switching on the approximation, we can use a smaller value of Zmax, since high multipoles 
will not have time to grow. 

The same treatment can be applied to ultra-relativistic species, which are identical to 
photons in this regime, except for the fact that they remain collisionless during reionisation. 
When the RSA is turned on for photons, it is better to follow ultra-relativistic species in the 
same way as photons, rather than with the fluid formalism described in Sec. |[ The RSA 
then removes three more differential equations, and cures the fact that the UFA turns out 
to be inaccurate at late time. Hence, the default version of CLASS treats ultra-relativistic 
species first with exact equations, then with the UFA (inside the Hubble radius and until 
photon decoupling), and finally with the RSA (inside the Hubble radius and after photon 
decoupling) . 

An expression for the RSA was discussed in the Newtonian gauge by Doran [jll| . Soon 
after, a somewhat simpler RSA (neglecting reionisation) was also introduced in the CAMB 
code, which uses the synchronous gauge. Here, we will derive an approximation comparable 



to that of Doran [11|, but valid in the synchronous gauge. 
4.1 Relativistic relics (massless neutrinos) 

We start with the simplest case, that of ultra-relativistic species ur. We combine the first 
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two equations of (|3.l|) into 



6r 



fr2 

" + —5 

r + 3 Our 



3 3 



(4.1) 



Inside the Hubble scale (i.e. when kr ^ 1) we can assume in first approximation that 
I Cur I ^ \^\it\ a.nd neglect the shear in the RSA. Also, since we are looking for a smooth 
(non-oscillatory) particular solution of this inhomogeneous equation, we can assume that 



<^ k \5ut\- We conclude that the RSA for 5ut is simply 



(4.2) 



Note that in the synchronous gauge, h' coincides with —26'^^, where Scdm is the cold dark 
matter density contrast. Deep inside the matter-dominated regime, Scdm oc a oc r^, so /i' is 
linear in r, and h" is a constant. The RSA for S^r is therefore nearly static. Concerning ^ur, 
its value in the RSA is given by the exact energy-conservation equation 5'^ 
0, namely 



3 Cur 



1 



h' . 



(4.3) 



In practice, we must extract h' and h" from the Einstein equations in the synchronous 
gauge, that read 



2k^rj - -h' 
a 



h" + 2-h' 

a 



2k^i]' 
- 2k^r] 



87rC/a^(5ptot, 



(4.4a) 
(4.4b) 
(4.4c) 



Here we do not need the fourth equation sourced by the shear. The difficulty comes from 
the fact that in order to infer (5uj- we should compute h" , and for doing that we need to 
combine the first and third equation, i.e. we need to know Sptot, which depends itself on 
(5ur- Fortunately, we can notice that if we omit 5ur in the computation of 5ptot, we make 
a tiny error, since during matter domination |(^/9ur| ^ I'^Pcdml- Hence, it is good enough 
to evaluate the first Einstein equation with 5ur set to zero. The same is not true for the 
second equation, since the synchronous gauge is comoving with cdm, so one has ^cdm = 
by construction. As a result, neglecting 6^^ in the computation of ^tot and rj' leads to a 
significant inaccuracy in the solution for r]. Hence, we introduce the following scheme: 



1. we compute Jptot assuming 5^ 
equation. 



0, and obtain 2k'^ri — —h' from the first Einstein 



2. using the fact that during matter domination |(5ptot| ^ |<5ptot|) we notice that 



h" + 2-h' -2k'ri 
a 



and hence to very good approximation 



-h' -2k'^7] 



(4.5) 



h" = -2-h' + 2fc2^ . 
a 



(4.6) 
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We then infer the fohowing RSA for Jur from Eq. ( [4 .21) : 



Snr = ^(-h'-k'v] ■ (4.7) 



a 



This formula is practical since r/ is one of the variables that we integrate over time, 
and h' has been inferred in the previous step. 

3. We impose the free-streaming solution for O^t (Eq. ([4.3D) and set a^T = 0. 

4. We use the remaining Einstein, continuity and Euler equations to evolve the system. 

4.2 Photons 

For photons, the solution is a bit more complicated since the baryon-photon coupling 
cannot be neglected during reionisation. We then need to find the particular non-oscillatory 
solution of (cf. Q) 

5% + = -\h" + ^A:^ -^{Ob- 0,). (4.8) 
6 6 6 6Tc 

Once again we will neglect the shear and search for a particular solution slowly varying 
with time (|(5"^| ^ |/c^5^|). In order to deal with the coupling term, we expand the solution 
in powers of r^^. The zeroth-order solution is exactly similar to that of massless neutrinos, 
Eqs. (^), ([4. 31) . The first-order solution should satisfy 

y'^7 = -^/^"-^(^^ + ^/^'), (4.9) 



in which h" can be replaced using Eq. ( |4.6| ). This approximation turns out to work out 
very well, unlike the zeroth-order solution. The velocity is then given by the exact energy- 
conservation equation 

= -\h' - . (4.10) 

We take the derivative of the previous result for 5-y, assuming that h" is time-independent, 
and using once more Eq. ([4.6|). We obtain: 



Tr. 



h + Ih'] + {eu Ih" 



(4.11) 



in which h" can be replaced using Eq. (|4.6|) . However the exact expression of 0^ depends 
again on 6^. Like before, we use a perturbative scheme in and replace 6'^ above by its 
expression at first-order in r~^. 
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4.3 Summary of RSA equations 

In summary, the RSA consists in neglecting 5^ and 5ur in the evolution of 5ptot in the first 
Einstein equation, which allows us to obtain h' , and then in imposing 



5^ 












l5ur 






- 4"' ■ 


<7ur 


= . 



+ 



(4.12a) 

+ 1 0b + - -/i' + fc^r? ) , (4.12b) 

(4.12c) 
(4.12d) 

(4.12e) 
(4.12f) 

This scheme is set to be the default one in CLASS, as long as the precision variable ra- 
diation_streaniing_approximation remains set to rsa_MD_with_reio. For comparison, 
some cruder schemes can be used: if the same variable is set to rsa_MD, the code will use 
the above expressions at zero order in (i.e, with 5^ = and 0^ = Q^-r). If it is set to 
rsa_none, the radiation perturbations are just set to zero. 

4.4 Comparison at the level of perturbations 

In Fig. ^, we show the evolution of ^-y, Q\i and r] between recombination and present 
time, for a particular wavenumber = 0.1 Mpc^^. We compare two RSA schemes with 
the exact evolution obtained by integrating all multipoles at all times. We always keep the 
UFA approximation of Sec. ^ turned off. We see that immediately after switching on the 
RSA, our approximation for (and also for ^ur, which is not shown) matches accurately 
the exact evolution averaged over a few oscillations. This would not be the case with 
several simpler RSA schemes which assume full matter domination and an exact linear 
growth of /i'(t). In the models used for the figures, reionisation takes place at z^^ = 10 and 
T = 4458 Mpc. It induces a clear feature in 5^ and Oh (having impact on rj) which is well 
captured by the terms proportional to in the full rsa_MD_uith_reio scheme. 

4.5 Comparison at the level of temperature/polarisation spectra 

The precision parameters governing the evolution of perturbations in the late universe are: 

• two parameters defining the time at which the RSA is switched on. For each wavenum- 
ber, we stop evolving photon and ur perturbations when the two conditions 



kr = r/Tk > radiatioii_streaming_trigger_tau_over_tau_k 



and 



Tc/t > radiation_streaming_trigger_tau_c_over_tau 

are satisfied (i.e., photons are sufficiently decoupled, and the wavelength is sufficiently 
deep inside the sub-Hubble regime). In principle, it would be possible to switch on 
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Figure 8: Evolution of the quantities 5^ (top left), 9-^ (top right), 6h (bottom left), 77 (bottom 
right) for the mode k = 0.1 Mpc~^, as a function of conformal time (in Mpc), between a time 
chosen soon after photon decoupling (or slightly before reionisation in the 6'b plot) and today. The 
red curves show the result of an exact integration with no truncation or approximation. In the blue 
curves, the RSA is turned on around r] = 470 Mpc, just before the solutions become unphysical 
due to the Boltzmann equation truncation at l_max_g=12, l_max_pol_g— 12, l_max_ur=28. For 
comparison, in green, we show the result for 5-^ when the terms in r^T^ are neglected in the RSA, and 
those for 0f, and 77 when the radiation multipoles are all set to zero instead of using a free-streaming 
solution. 



the RSA at different times for photons and ur species, but for simplicity we did not 
consider this option. 

• l_max_g, l_max_pol_g and l_max_ur define the number of photon temperature, 
photon polarisation and ur multipoles which are integrated until the RSA is switched 
on. 

In order to compute some reference spectra to be used throughout this section, we 
fix the precision parameters according to the file cl_permille .pre, increase l_max_g, 
l_max_pol_g and l_max_ur to 3000, and choose such large values of the trigger parameters 
that the UFA and RSA are never employed. We wish to compare these reference spectra 
with those from runs with/without the RSA, in which l_max_g, l_max_pol_g and l_max_ur 
are kept fixed to reasonable values. We choose l_max_g and l_max_pol_g to be equal to 
18. Since we do not want to use the UFA approximation in this comparison (in order to 
focus only on the impact of the RSA), we fix l_max_ur to a larger value, namely 50. This 
setting, called no-rsa in Table ^, leads to a 0.01% error both in the temperature multipoles 
and in the matter power spectrum for k < l/iMpc~^, as shown in Fig. ^. Finally, in the 
run called rsa, we switch on the default RSA scheme, with the trigger values specified in 
Table |^. With such settings, the error in the temperature (and also polarisation) multipoles 
remains as small as without the RSA, while the error in the matter power spectrum grows 
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moderately to 0.04% (see Fig. P). However, the running time is reduced considerably, as 
shown in Table |5|: with the Runge-Kutta integrator, the RSA leads to a 66% speed up. Note 
that the ndf 15 remains much better than the Runge-Kutta integrator when the RSA is 
employed (by a factor 2) while it experiences difficulties in following oscillatory solutions in 
absence of a RSA. However, the combination of our RSA scheme and ndf 15 integrator leads 
to very nice performances (speed-up by a factor 4 with respect to Runge-Kutta without 
any RSA). 



precision setting: 


reference 


no-rsa 


rsa 


l_max_g 


3000 


18 


18 


l_max_pol_g 


3000 


18 


18 


l_max_ur 


3000 


50 


50 


radiation_streaming_trigger_tau_over_tau_k 


oo 


CXD 


100 


radiation_streaming_trigger_tau_c_over_tau 


oo 


CXD 


2 


ur_f luid_trigger_tau_over_tau_k 


oo 


OO 


oo 



Table 4: Three settings for the parameters governing the Boltzniann truncation and RSA. All 
other parameters are fixed with the file cl_permille .pre of the public CLASS distribution: in 
particular, radiatioii_strecmiing_approximation is set to rsa_MD_with_reio. The reference run 
never uses the UFA and RSA, and cannot be affected by the Boltzmann truncation. The second 
and third settings share the same truncation multipoles, and differ only by using the RSA or not. 




Figure 9: Temperature and matter spectra for the runs rsa and no-rsa, normalized by reference 
spectra. Corresponding precision parameter settings are described in Table ^. Using the RSA with 
such settings leads to equally accurate CMB spectra, and to a very small error on the matter power 
spectrum. The running time is however reduced considerably. 

In CAMB, the RSA is somewhat cruder, since it neglects reionisation and uses an explicit 
cosmology-dependent relation giving h' an h" in terms of pi,, pcdm, ^h-, <^cdin and /c, valid only 
deep in the matter-dominated regime. We present in a companion paper the comparison 
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Runs with/without a RSA: 


no-rsa 


rsa 


rk 

ndf 15 


55.8s 
84.4s 


33.6s 
14.5s 



Table 5: Execution time of the perturbation module in seconds, with the precision parameters of 
the file cl_permille .pre, modified as described in Table |[ 

between the matter power spectrum P{k) computed by CAMB and CLASS. In order to get 
an accurate P{k) with CAMB, one is forced to deactivate the RSA approximation (called 
late_rad_truncation), precisely for the above reasons. Our scheme does not lead to 
a significant error on the P{k), and is more model-independent: it involves only metric 
perturbations and does not need to be modified in the presence of other components playing 
a role during matter domination (e.g. with warm dark matter or early dark energy). 

5. Conclusions 

In Fig. we summarise the different approximations used by the code in the {k, r) plane 
when computing the C/'s up to 3000 and the P{k) up to l/iMpc~^. The figure corresponds 
to the precision settings of the file cl_3perniille .pre. Exact equations are used only in 
the band corresponding to Hubble crossing for each mode, as well as in the super-Hubble 
region with non-tightly-coupled photons, in which all quantities evolve very slowly and 
integration is very fast. 

With these approximations, for ACDM, the perturbation module only spends a sig- 
nificant time in the region corresponding to Hubble crossing for each mode, and to the 
sub-Hubble evolution before photon decoupling. At early times, stiff equations are avoided 
thanks to the Tight-Coupling Approximation. Well-inside the Hubble radius and until pho- 
ton decoupling, the different modes oscillate and integration is time-consuming: however, 
the number of equations is kept small (of the order of 30 in total) thanks to the UFA. After 
photon decoupling, the code only needs to integrate over 4 equations with very smooth 
solutions, and the time spent by the code in the RSA region is negligible. 

Other approximation schemes can be introduced in more general cosmological models. 
The case of massive neutrinos and non-cold dark matter relics is discussed in a companion 
paper. More exotic cases may require further approximations which can be introduced and 
discussed case-by-case (CLASS is coded in such way that introducing a new approximation 
is as structured, codified and simple as introducing new species Q). However, the fact 
that CLASS uses an original stiff integrator means that for several purposes (as for the 
generalisation of TCA), new approximation schemes are not even strictly necessary. 
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Figure 10: Summary of regions in (k, r) space where ttie various approximations are used. The 
precision settings are taken from the cl_3permille . pre precision file which ensures a 0.3% precision 
on the Ci's till I = 3000. The full set of exact equations is used only in the white region. 
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A. Stiff integrator 

The standard numerical method for solving Ordinary Differential Equations (ODEs) is to 
use an adaptive step size Runge-Kutta solver. While this method is fast and accurate 
in simple cases, it may fail completely (or take a very large number of steps) when the 
problem is stiff. Stiffness occurs when at least two times scales of evolution in the problem 
differ substantially. A well known example is the Boltzmann equation in cosmology. If a 
distribution is kept in equilibrium by the (rapid) interaction with a background species, 
and we are interested in the evolution of the distribution on cosmological time-scales, a 
Runge-Kutta solver will oscillate around the equilibrium solution by means of very small 
time steps related to the time-scale of the interaction. 

This problem exists, for instance, in the early universe where baryons and photons 
are strongly coupled. One solution is to substitute the equilibrium solution into the equa- 
tions and thereby reducing the system of equations and removing the stiffness: this is 
the basic principle of the Tight Coupling Approximation (TCA) discussed in Sec. H This 
approximation removes the stiffness of the Boltzmann equation in the vanilla scenario. 

However, stiffness may easily be reintroduced by trying to incorporate new physics 
into the code. The proper TCA approximation should then be derived and implemented 
in the code, which would require a great deal of familiarity with the code from the user. 
In CLASS, however, the user can just take advantage of the implemented stiff solver. 

The ndf 15 algorithm is a variable order (1-5) adaptive step size solver based on the 
Numerical Differentiation Formulas of order 1 to 5. The step size is adaptive but quasi- 
constant, meaning that the formulas used are based on a fixed step size. Each time the 
step size changes, the code will update the backward differences by interpolation to reflect 



this new step size. The algorithm is described in [17| 
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Whereas Runge-Kutta methods are exphcit, meaning that the next step can be com- 
puted directly from the previous step by elementary operations, ndf 15 is a fully implicit 
method. This means that at each time step, we must solve a system of algebraic, possibly 
non-linear equations, which is accomplished by Newton iteration. This requires a numer- 
ical computation of the Jacobian and the solution of systems of linear equations. As it 
is evident, all this can make each step very expensive, so a number of strategies must be 
implemented: 

• Reusing Jacobians: 

The Jacobian usually changes more slowly than the solution itself, so an attempt is 
made to reuse the Jacobian - it will only be recomputed if Newton iteration is too 
slow. 

• LU decomposition: 

The same linear system must be solved repeatedly with different right hand sides, so 
we should of course store an LU decomposition. Because the system of equations can 
be large, O(IOO), we need sparse matrix methods for this. 

• Backward Interpolation: 

Since the method stores a matrix of backward differences, it is fast to infer the value 
and the derivatives of the solution at points before the current point by interpolation. 
We only need the values of some of the components to calculate the source functions, 
so we do not need to interpolate the rest of them. 

When the number of equations is larger than about 10, it is advantageous to use sparse 
matrix methods, and if the system is somewhat larger, the difference in execution time can 
differ by orders of magnitude. Usually sparse matrix methods are developed in order to 
save both time and memory, but for our purposes only execution time matters. We created 
a small sparse matrix package for our purpose based on jl^. Some important features of 
this package are: 

• Column Pre-ordering: 

The matrices which appear are of the general form I — cJ where / is the identity 
matrix, J is the Jacobian and c is some constant. Since the Jacobian is close to being 
structurally symmetric (if yi couples to yj, it will often be the case that yj couples 
to yi), it is advantageous to use the Approximate Minimum Degree column ordering 
of the matrix C = + J. Using this pre-ordering reduces the number of non-zero 
elements in the corresponding L and U factors by a factor of a few, which leads to a 
sizable reduction in the time needed to solve the linear systems. 

• LU re-factorisation: 

If we have already calculated a LU factorization for some Jacobian, which is struc- 
turally identical to the current Jacobian, we can use information saved during the 
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factorisation of the former to factorise the new matrix in a fraction of the time. 
Specifically, we store the pivot ordering and the reach of all the sparse right hand 
sides used in forming the LU-decomposition. 

• Fast Jacobian Calculation: 

If the same pattern of the Jacobian is found repeatedly, we can use this pattern to 
speed up the calculation. Taking advantage of the sparsity of the Jacobian, we can 
group the columns together and form the Jacobian using only a fraction of the usual 
n function evaluations, n being the number of equations. 

The CLASS user can choose to use the Runge Kutta or ndf 15 integrator by switching 
the precision parameter evolver to either rk (=0) or ndf 15 (=1, default setting). To 
illustrate the power of ndf 15 in stiff situations, ndf 15 was less than 10% slower when the 
TCA was turned off as early as in the reference run used in Fig. |2[ As a comparison, the 
standard Runge-Kutta integrator was more than 10000 times slower than ndf 15 for the 
same model. This particular run represents an extreme case, but throughout this work we 
have presented various examples in which the ndf 15 performances are very good. 



B. Derivation of fluid equations for ultra-relativistic relics 

The goal of this Appendix is to establish the validity of the approximate shear derivative 



equation (|3.2D , which allows to treat collisionless species as an imperfect fluid governed 
by may less equations. In the future, our results could be used for computing higher 
order terms in (^]^), or more generally for understanding various properties of the linear 



perturbations of ultra-relativistic species. A discussion similar in spirit was presented in |19| 
for massive neutrinos, although the goal of that paper was to introduce a sharp truncation 
at / = 3, while we are searching for a truncation scheme that would take into account the 
transfer of power to higher Vs. 

B.l Formal solution 

Sticking to the notations of Ma & Bertschinger, the perturbations of ultra-relativstic species 
is described by a function -F(k, ^u, r) obeying to the collisionless Boltzmann equation 

drF{k, /X, r) + ikfiF{k, /i, r) = S{k, /x, r) , (B.l) 

where S stands for the gravitational source terms. The most general formal solution can 
be written as 



F(k, II, r) = F°(k, ii)e-'''^'^ + / e''''^'^^-^^ S{k, f)df . 

Jo 



(B.2) 



The initial function -F''(k, /j,) depends on the considered gauge and type of initial conditions. 
For instance, in the synchronous gauge, = for the growing ADiabatic (AD), Baryon 
Isocurvature (BI) and Cold Dark Isocurvature (CDI) modes; F^ has a non-zero monopole 
term for Neutrino Isocurvature Density (NID) initial conditions; and a non-zero dipole term 
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for Neutrino isocurvature Velocity (NIV). These statements can be checked from ref. |20|. 
Hence, in all cases, we can write 

FO(k,/i) = C7NiD(k) - ikfiCmvO^) ■ (B.3) 

More fundamentally, the fact that only has monopole and dipole contributions can be 
justified by the fact that neutrinos were initially in thermal equilibrium, forming a fluid 
with no anisotropic pressure or higher momenta. By causality, this remains true at any 
time after decoupling on super-Hubble scales. So, we are sure that the above form of F° 
is completely universal. In the synchronous gauge, the source reads 

S = -^h'-^{h' + 67j')P2{p) . (B.4) 

Thanks to a few integration by part, it is possible to absorb the // dependence, in order to 
be able to write this solution in Legendre space. The result is 

F(k,/i,T) = F°{k,fi)e-'^''^ (B.5) 

+^ e'^''^^-^^ {2k^v'ir) + h"'{f) + 67?'"(f )} df (B.6) 

-^{h" + 67]" -ikfi{h' + 6rj')} (B.7) 

+^e-"'^'^ {h" + 6i]" -tkfiih' + 6r,')}^^^ . (B.8) 

The last bracket contains the initial value of (/i' + 6r/') and (h" + 6r]"). The former vanishes 
for all types of initial conditions excepted NIV; the latter is non-zero for AD, NID and 
NIV. All initial condition terms in the first and last lines can be grouped and represented 
by two coefficients a and /?: 

Q(k) - ifi/3ik) = FO(k, ^) + ^{h" + 67]" - ikii{h' + 677')},=o (^-9) 
= l^ur + ^{h" + 6r/")} ^ - + ^(^' + 6^)} ^ , (B.IO) 

with e.g. (a, /3) = (20/(15+4i?ur)) 0) for adiabatic initial conditions (as can be checked from 
pO[] ), i?ur being the fractional contribution of ultra-relativistic species to the background 
density. It is now easy to expand the full solution in Legendre coefficients, using the 
definition F[^) = + l)FiPi{ij) and the fact that plane waves can be expanded 

in terms of spherical Bessel functions: 

FKk,r) = a{k)ji{kT) + mfiikr) (B.ll) 

+^2 (^(^ - ^)) {2fc'V(T) + h"'{f) + 6rj"'{f)} df (B.12) 

-^[ih" + Qv")Sm + l{h' + 6v')Sii^ . (B.13) 



The terms in the ffist line show how initial conditions propagate to later times, by just 
free-streaming. The other terms show how perturbations adjust themselves to the power 
injected in the system at any time by metric perturbations. 
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B.2 Sub-Hubble approximation 

Well inside the Hubble radius, the above formal solution can be approximated by a simpler 
expression. The results of this subsection are never used in our UFA scheme or in CLASS, 
but we present them for completeness, and also because the approximation performed in 
the next subsection will follow the same logic. 

The second line of the solution contains a convolution between a function which varies 
smoothly over a Hubble time (at least for kr ^ 1), and a Bessel function ji{x) with 
X = k{T — f) which oscillates over = 1/k. Bessel functions ji{x) peak near Xpeak = ^ + 1/2 
(in fact this statement is accurate only for very large I; for instance, ji{x) peaks near 
^^peak = 2 and j2{x) near Xpeak = 3.5). The integral on f G [0, r] corresponds to x G [0, kr]. 
The goal of this subsection is to find an approximation for this convolution for small / 
values. 

As long as /cr < 1, it is difficult to find a low-Z approximation for the convolution; 
the result is a function oscillating over a characteristic time = 1/k. In this regime, the 
integral brings an extra oscillatory contribution to the term a(k)ji{kT) + /3(k)j[(A;r); this 
explains while around Hubble crossing, the numerical solution for i*)(k, r) exhibits irregular 
oscillatory patterns, with very different peak amplitude between two consecutive periods. 

However, when kr ^ 1, the integral runs over a large range x G [0, kr]. For low /, this 
means that the convolution picks up significant contributions only near x = Xp^ak ^ kr. 
Near this value, the slowly-varying argument can be approximated as a constant, to be 
evaluated around f = {kr — Xp^akj/k, i.e. in very good approximation near r, since kr » 
a^peak- So, we can write: 



jl {k{T - f )) {2er^'{f) + h"'{f) + Qr,"'{f)] df 

{2k'^r]'{T) + h"'{T) + 6?7'"(r)} [ ji {k{T - f)) df . (B.14) 

Finally, still in this limit A;r ^ 1, the last integral can be approximated by 

1 /"^ . X yiT(//2 + l/2) 

- / jl ix) dx = , \i, \ B.15 
k Jo 2kT{l/2 + l) ^ ' 

A more accurate approximation scheme would lead to extra contributions involving time 
derivatives of the quantity between brackets in ( B.14| ): for instance, the next order term 
would be of the type 

[2k\"{T) + /iW(t) + 6r?(4)(T)} ^ (B.16) 

with 7 being a coefficient of order one. However, since inside the Hubble scale metric 
perturbations vary over a Hubble time r ^ r^, we can keep only the leading source terms 
with the highest power of k: 

aji{kr)+(3jl{kr) + ^^v' , (B.17) 

ajiikT) + (3j'i{kT) - ^h' , (B.18) 
aji{kT)+Pjl{kT) + ^rj' . (B.19) 











3k^^' 




2<7m 


= F2 
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In the code, we do not use directly these asymptotic approximations, and try instead to 
close the system of differential equations with a trunction at order two, like for a viscous 
fluid. 



B.3 Exact truncation formula 

First, in order to manipulate more compact equations, we define: 

Giik,T)^Fiik,r)-^!^{h" + Qr]")5io + ^{h' + 6r^')6n^ . (B.20) 

Since j'i{x) = ji-i{x) — ^-^ji{x) and Gi contains a term jiikr), we suspect that it obeys 
approximately to a similar relation. We use equation ( B.13|) for computing exactly the 
difference G'l — kGi-i + ^^Gi, which would vanish if only the term aji{kT) was contributing. 
We use the identities 

j'l'ix) = jti(x) - ^-^j'lix) + ^-^Mx) (B.21) 

and ji{0) = 6io. We find 

G'l - kGi^i + = P^^3i{kT) + ^ {2k^7i'{T) + h"'{T) + 6r?'"(T)} 6io 

2 



r K{t, f ) {2k'^7]'{f) + h"'{f) + 6??"'(f )} df (B.22) 
Jo 



where we defined 



K{t, f ) = kjl {k(T - f )) - kjl^i (k{T - f )) + {k{T - f)) (B.23) 



This expression simplifies to 



K{t, f) = J-±l[—^\ji {k{T - f)) (B.24) 



T — T 



This trunction scheme is not pratical in general. However, the goal of the UFA is to find a 
way to close the system at low / not at all times, but only deep inside the Hubble radius. 

B.4 Sub-Hubble truncation formula 

In the limit /cr ^ 1, equation ( [B.22| ) can be simplified for two reasons. First, the Bessel 
function varies over a time scale = l/A; ^ r, so 

j'likr) » ^ (B.25) 

Hence, in this limit, the last term in the identity ( p^ ) can be omitted, which implies that 
the first term (proportional to /3) in eq. (B.22) is always negligible, even in the presence 
of isocurvature modes. Second, we can devise an approximation for the convolution, by 
noticing once more that it involves metric perturbations wich vary smoothly over a Hubble 
time (at least for kr ^ 1), and the quantity ^^~^ ji{x) with x = k{T — f) which oscillates 
over a period of order = l/k. 
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When kr ^ 1, the integral runs over a large range x G [0,A;r]. For low /, this means 
that the convolution picks up significant contributions only near x = Xpeak ^ kr, while in 
this range ji{x) ~ ^ji{x). Near f = {kr — Xpeak)/^ — ''") the slowly- varying metric 
perturbations can be treated as a constant term. So, we can write: 

jl {k{T - f )) {2k'r^'{f) + h"'{f) + 6??'"(f )} df 



T 



{2erj'{T) + h"'{T) + 67?'" (r)} (737) Jl (Hr - r)) df (B.26) 
(737) (B.27) 



with 



Let us estimate the error made in these approximations. Since we are inside the Hubble 
radius with smoothly varying metric perturbations, and since the integral in eq. ( [B.26| ) is 
of order r, the leading term in eq. ( [B.26| ) is of order {Jt^Trj'). The other terms of order 
(Tr]'") and {rh'") can be neglected. If instead of considering r]'{f) as a constant we perform 
a Taylor expansion of this function around r, we find that the next order contribution to 
eq. (B.26) is of order (krr]"). Finally, an explicit calculation shows that the approximation 
performed in ( [B.27| ) amounts in neglecting terms of order k~^ with respect to terms of 
order r. In summary, we obtain the following approximate truncation equation: 



^ + 1^ 2 



G'l - kGi.i + —^Gi = {2k'r,'{T) + h"'{T) + 6??'"(t)} 61 







_4(, + l),/f + . (B.28) 



For / = 2, the integral is equal to -1/3, since {ji{x)/x)' = j2{x)/x and lim^_>.o[ji(x)/x] = 
So, 

G'2-kGi + ^G2 = ^v' + o{^,^ . (B.29) 
Replacing the Gis by the appropriate momenta, we find: 

< = -\'^ur + \0^r - \{h' + 677') + 2ri' + ^ . (B.30) 

The last term 2r]' comes from our approximation for the convolution. Without a full 
treatement like the one presented here, one would miss this term and obtain the truncation 
formula called uf a_mb in the code, based on simply assuming Gi oc jiikr). However, with 
this extra contribution, the two terms in 77' cancel each other, and at leading order we end 
up with 

a'^, = --a^, + leur-lh' , (B.31) 
T 6 6 

which is precisely what we call the uf a_class truncation scheme in CLASS. 
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